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ABSTRACT 



We use path integrals in order to estimate merger rates of dark matter haloes 
using the Extended Press-Schechter approximation (EPS) for the Spherical Col- 
lapse (SC) and the Ellipsoidal Collapse (EC) models. 

Merger rates have been calculated for masses in the range lO lo M h^ 1 to 
lO 14 M h _1 and for redshifts z in the range to 3. A detailed comparison between 
these models is presented. Path approach gives a better agreement with the ex- 
act solutions for constrained distributions than the approach of |Sheth fc Tormen 



| (2002). Although this improvement seems not to be very large, our results show 
that the path approach is a step to the right direction. Differences between the 
two widely used barriers, spherical and ellipsoidal, depend crucially on the mass 
of the descendant halo. These differences become larger for decreasing mass of 
the descendant halo. 

The use of additional terms in the expansion used in the path approach, other 
improvements as well as detailed comparisons with the predictions of N-body 
simulations, that could improve our understanding about the important issue of 
structure formation, are under study. 



Subject headings: galaxies: halos - formation; methods: analytical; cosmology: large 
structure of Universe 
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Introduction 



The development of analytical or semi-numerical methods for the problem of structure 
formation in the universe helps to improve our understanding of important physical 
processes. A class of such methods is based on the ideas of Press & Schechter ( 1974[ ) and 
on their extensions. These extensions are called Extended Press-Schechter Methods (EPS), 



and they are presented in the pioneered works of Peacock & Heavens (1990), Bond et al. 



(1991) and Lacey & Cole (1993) 



In this section we summarize useful relations about density perturbations, filters, stochastic 

processes and path integrals that are used in the following calculations. 

The overdensity at a given point r of the initial Universe is given by the relation: 

p(r) - p b {r) 



8(v) 



(1) 



In the above relation, p(r) is the density at point r of the initial Universe. Index b denotes 
the density of the background model of the Universe. The smoothed density perturbation 
at r is defined by the relation: 



5(r;R) = j5{x)W(r-x; #)d 3 x 



(2) 



where W is a filter with characteristic radius R. Using the convolution theorem to transform 
both sides, we have: 

5{k;R) = 5{k)W{k;R) (3) 



In our calculations we use the sharp in k space filter given by: 

1 3(sinx — xcosx 



W KS {r-R) 



6n 2 R s 



X" 



x = r/R 



W KS (k;R) = H(l - kR) 



(4) 
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where x = r/R and i7 is the Heaviside step function denned by: 



H(x) 



0, x < 
§, x = 



(5) 



1, x > 

For spherically symmetric kernels, such the ones we examine here, the variance of the 
overdensity at scale R is given by : 



S(k f ) = o*(k f ) 



A 2 (k)\W(k;k f )\ 2 d\nk 



(6) 



where kf = 1/R and A 2 = j^gPs(k)- P s {k) is the power spectrum . For the k-sharp filter 
the evolution of smoothed 5 as a function of S is governed by a Langevin equation: 



d5 
dS 



n(S) 

<n{S)n{S') >=5 D {S-S') 



(7) 



(Langevin (1908), Coffey et al. (2004)), where 5d is the Dirac delta function, that 
describes the following Markovian stochastic process: Trajectories start from the same 
position (So, S ) at the (S, S) plane and evolve according to Eq.7 as S decreases. S is 
completely analogous to time t in ordinary problems involving stochastic processes. The 
probability P(S,8/So,So)d6, a trajectory that starts from the position (So, So) passes at S 
from a value of 5 in the interval [5, 5 + dd] satisfies a Fokker - Planck equation (Coffey et al 
2005), that leads to the diffusion equation: 

dP(S,5/S ,5 ) Id 2 



OS 



2d5 2 



[P(S,S/So,So)} 



with solution: 



P(S,S/So,S 



: exp 



(5 - So) 2 
'2(S-S ) 



(9) 



(10) 



'2ir(S - So) 

In Sect. 2 we give the basic equations from the path integral approach method. In Sect. 3 
first crossing distributions and structure formation are discussed. In Sect. 4 analytical 
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relations for merger rates and mean merger rates are presented. Results for different models 
as well as comparisons are presented. In Sect. 5 we summarize the method, the results and 
we give a short discussion. 



2. Path integral approach: Basic Equations 



Path integrals are power tools for the study of various fields of theoretical physics, and 
they have been studied extensively in the literature, (Wiener ( |1921 ) Feynman & Hibbs 



(1965), Grosche & Steiner (1998), Chaichian & Demichev (2001)). They are given by the 
sum over all possible paths, satisfying some boundary conditions. Applications appeared 



for the problem of structure formation in the Universe from Bond et al. (1991), Maggiore 



& Riotto (2010) and Simone et al. (2011). We follow the formalism of the above authors. 



According to the above described picture for the formation of structures, we consider an 
ensemble of trajectories all starting from the point (So, So) and we follow their evolution for 
the "time" interval [So, S]. The interval is discretized in steps AS such as Sk = So + ke where 
k = 1,2. ..n and S n — S.A trajectory is defined by the points (So, So), (Si,Si)...(S n ,S n ). The 
probability density the variables S(Si),S(S2)-.-S(S n ) take the values 5i, <5 2 ...<5 n respectively 
is: 

W(So,S ;S n ,S 1 ,S 2 ...S n ) = (S^S^ - S 1 )...S D (S(S n ) - S n )) (11) 



Using the integral representation of the Dirac delta function : 

^+oo dA 



S D (x) 



-iXx 



2tt 



(12) 



and the fact that each one of the variables S(Si),S(S 2 )...S(S n ) follows a central Gaussian 
distribution, the probability Il e of arriving at (S n , S n ) starting from (5*0, S ) using discrete 
trajectories of step e that never exceeded a barrier B that depends on S is: 



n e (S'o, So] S n , S n ) 



dS[ 



dS' 2 ... 



Br, 



dS' n _ 1 W(S , S ; S n , S[, S' 2 ...S' n ) 



(13) 
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where, 



and 



W(S ,5 ; S n ,5' 1 ,5' 2 ...5' n 



+oo dAi r+oo^ f+°° d \ n 



-oo 27r J_ 00 2tt J_ 00 2-n 



(14) 



n n ^ n n 

z' = i/ J h(B k - B n ) ^k5' k - - 2j E ^jSkj 

k=l k=l k=l j=l 

Expanding in a Taylor series the quantity e l ^=i x k(B k -B n ) ag we jj ag ^ e difference 



(15) 



Bk — B n around S n (see Lam & Sheth (2009) and the discussion therein), we have 
W = W<® + + + where 



W {0) 



—— / V\ exp[i hh ~ g E E AfeA A 

^ ' J ~°° fc=l fc=l 7 = 1 



(16) 



(2tt)' 



E VXX k exp[i ^ hh - 2 E E 

fc=i ^~ 00 fc=i fc=i 7=1 



(17) 



n 11 r+oo n 1 n n 

W {2) = E E ^ / DXXkXj exp[z £ A*fc " 2 E E A * A ^ 

\ ' k=l j=l 00 fc=l fe=l j=l 



where we set: 



and: 



+00 r+oo r+oo 

VX= dAi.. / dA„ 

00 J —00 J— 00 

+00 ^(p) 

Cfc = / . — (Sk — S n ) p 
z — ' p 

p=l 



+ OO +OO (p) (g) 

p=l <,=! ^ 



The symbol denotes ^§®/s=s n . 
It is straightforward to show that and are connected to W^°\by: 



(18) 

(19) 
(20) 

(21) 



(22) 



k=l 



k=l j=l 
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where = ^ and d^d^.-d^ = gSk ® k — ^7. The quantity 1/0°) is known as Wiener measure 
and is given by: 

1 1 - 

WW(5 , 5 = 0; S n , 5 1 ,6 2 ...6 n ) = — — - 2 exp[-- ]T(4 - 4-i) 2 ] (23) 

^ 7re ^ e k=i 

We note that the integral f^ n ^ U t= o(S ,5 ; S n ,5 n )d5 n equals to the number density of 
trajectories that have never crossed the barrier. This number density is a decreasing value of 
"time" S, since the number of trajectories that pass the boundary increases with increasing 
S. However, the rate of change of the above integral shows the number of trajectories that 
cross - for first time- the barrier at S. This is the first crossing distribution J 7 that satisfies 
the Eq. 

Q fB(S n ) 

F(S n ) = -Qg-J n e=0 (S , 5 ; S n , 5 n )d6 n (24) 

Assuming that every component (16), (17) and (18) satisfies a diffusion equation, we can 
write: 



F {j) (Sn) = -~Il { £ (So,5 ;S n ,5 n )\s n=Bi s nh J=OA 



(25) 



Finally we have: 



^ (0) (Sn 



B n — 5 



+00 



^ (1 »(«=^W.-«5o)E i #4 ri 

p=l p - 



27T(5 n - S0) 3 / 2 



exp 



exp[ 



[B n -5 
2(S' ri — So 



2 1 



s (Si-S )^ 2(S l -S 



dS t (26) 



that in terms of the confluent hypergeometric function U ( Abramowitz & Stegun (1970)) 
can be written: 

+00 



^\s n ) = i\M e -^ j2(s n - s y- 2 - 



2n 



lY B^Y (p-l)u(p- 



P =i 



1 3 

2' 2' 



~-j= - S r 3 ^B^T (p-1)u(p-1,\^) (27) 
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where ip = z$~) an d F is the complete gamma function. 
Thus T(S n ) = J* 0) (S n ) + F m (S n ) is written in the form: 

Tp ath (S, B(S)/S Q , 5 Q ) = ^ e { J_ So)V2 G ^ 5 > 5 o) ( 28 ) 

where the quantity G is defined by the relation: 

+00 



G(S, S , B, So) = B(S) -5 + E ^ )P B {p \S)T(p - \)U (p 



(29) 



In the above relation we have set S = S n , B(S) = B n , U(0,a,b) = 1 and we rewrote the 
arguments in a different form, useful in what follows. 

It is also useful for our purposes to consider a barrier that is a function of redshift 
z too, B = B(S,z). Additionally, we also assume that the value So satisfies the 
relation S = B(S ,z ). However, the constraint (S, B(S, z)/S , S ) is equivalent to 
(S, B(S, z)/S , B(S , z )) or without any confusion (S, z/S , z ). Thus, given that a 
trajectory crosses, for first time, the barrier B(So,z ) at redshift zq, ^(S, z, /So, zq) gives 
the probability this trajectory to cross, for first time, the barrier B(S,z) at z. We also 
denote ^(Sjz) = F{S,z/Sq — 0,z — z in ), where z in satisfies B(0,z in ) = 0. 
Obviously, in the above relations B(S) has to be replaced by B(S,z) and B^ P \S) by 
B^\S,z) = j^B(S, z). Thus, it is more convenient to write: 



G(S, So, z, z ) = B(S, z) - B(S , z ) + ^ ( ^' flM(g, z)T(p - l -)U (p - 1, (30) 

The use of the above formula in practice, requires taking account a finite number of terms. 
Due to the behavior of the confluent hypergeometric function, it is not obvious that the 
leading order term dominates the sum. In our calculations, we tried an increasing number of 
terms N T (thus the sum above extends from p — 1 to p — N T ). We found that for N T > 12 
the results are the same. So, the value of Nt = 12 is sufficient for our purposes. We note 
that we also checked the results for very large number of Nt as 100 and no differences were 
found. 
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3. First crossing distributions and structure formation 



First crossing distributions are connected to structure formation (Bond et al. (1991), 



Lacey Sz Cole (1993)). We consider the variable M that is the relative excess of mass at 



scale R. This is written in the form: 



M(R) — M b (R) 3 ( fi tM2l 
M ( R ) = = ™ / S ( r ) rdr (31) 



M b (R) R3 J 

where M(R) is the mass contained in a sphere of radius R of the Universe and 
M b (R) = |7rp feJ R 3 is the mass contained in a sphere of radius R of the unperturbed model of 
the Universe that has a constant density p b . 

It is easy to check that in the case 5 follows a central Gaussian, then M follows a central 
Gaussian with the same variance. Thus, it is reasonable to assume that the relation between 
S and R (see Eq. 6) can be transformed into a relation between S and M. This assumption 
should have a complete physical meaning if the volume associated with the filter was that 
of a sphere of radius R. The volume associated with a filter F is given by: 

4 r+°° 

V f = -tiR 3 47iW F (r;R)r 2 dr (32) 

The Gaussian filter has an infinite extent and so it is difficult to understand the physical 
connection between R and M and for the k-sharp filter the above integral does not exist 
at all and so the volume is not even well defined. However, since the k-sharp filter is so 
convenient for the analytical approximation of the problem studied here, we followed the 
usual procedure that is to assume a mass M = ^7rp b R 3 associated with the k-sharp filter. 



We note that in Lacey & Cole (1993) a volume V = 6n R is quoted for the k-space filter, a 



result that is used without much justification, see the details in Maggiore & Riotto (2010). 
The second point is to connect the first crossing distribution of trajectories with the number 
density of haloes. This is done by using the following argument: The probability a mass 
element at redshift z belongs to a halo of mass in the range M, M + dM denoted by 
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f(M, z)dM equals to the probability a trajectory crosses, for first time, the barrier B(S, z) 
between S, S + dS denoted by F(S, z) \ dS |. Variables S and M are connected by the 
relation S = a 2 (M). The described equation can be written in the form: 

f(M, z)dM = F(S, z) | ^ \dM (33) 

and since we assume that all trajectories start from the point (So = 0,5o = 0) this is an 
unconstrained probability For the constrained case, we write: 

f(M, z/M , z )dM = F(S, z/S , Zo)\j^\ dM (34) 

A form of the barrier that results to a mass function that is in good agreement with the 
results of N-body simulations is B E c(S, z) given by: 

B A°f^\ = 1 + r J, „ (35) 
^cBsc(z) [aB 2 sc (z)p 

In the above Eq. a, (3 and 7 are constants. The above barrier represents an ellipsoidal 

collapse model (EC), |Sheth, Mo fc Tormen | ( pOOlj ), |Sheth fc Tor"ni^n~| ( [2002| (ST02 



hereafter). The barrier depends on the mass (S = S(M)) and it is called a "moving 
barrier". The values of the parameters are a = 0.707, /3 = 0.485, 7 = 0.615 and are 
adopted either from the dynamics of ellipsoidal collapse or from fits to the results of N-body 
simulations. In the above relation B sc (z) = 1.686/ D(z), where D(z) is the growth factor 
derived by the linear theory, normalized to unity at the present epoch. B sc (z) is the barrier 
for the spherical collapse model that results from the above relation for a = 1 and (3 = 0. In 
the plane (S, S) the line 5 C = B S c(z) is parallel to the S axis. The physical picture is that 
in an Einstein-de Sitter Universe a spherical region collapses at z if the linear extrapolation 



of its initial value 5i n up to the present epoch equals to Bsc( z ) (see for example Peebles 



(1980)). 



The above relation can be written in a more convenient form: 

B EC (S,z)=p c (z) + q c (z)S^ (36) 
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where p c (z) = WiB sc (z),q c (z) = w 2 B^ 3 (z) with W\ = y/a,w 2 = (3a°^~ 1 ,w 3 = 2j — 1. For 
the spherical collapse (SC) model, T can be calculated analytically. The analytical solution 
for the SC model is given by: 



?sc{S, z) 



Pc{z) 



exp 



2S 



(37) 



Sheth & Tormen Q2002|) have shown that a good approximation of T for the EC model is 

B 2 EC (S,z) 



given by the relation: 

Fec-st^S, z) 
where 



V2ttS 3 



T(S, z) | exp 



2S 



T(S, z) = B EC (S, z) + J2 { -JT^-k B ^c(S, z 



fc=i 



We note that for the constrained case the above relations are written: 



J 7 sc(S, zj So, Zq) 



B 



SC 



z) - 5n 



v/27r(S-5 ) 3 



exp 



[B 



sc 



Z) - 5n 



2(S - S 



(38) 



(39) 



(40) 



and 

FeC-St{S, zj Sq, Zq) 

with 



1 



2tt(S - S ) 3 



T(S,z/S ,z ) | exp 



[B ec (S,z)-5q] 2 - ] 
2(S-S ) 



T(S, z/So,z ) 



B EC (S, z)-6o + J2 ^k B *c(S, z) 



(41) 



(42) 



Comparing the expressions (28) and (41) describing the path approach and the ST02 
approximation respectively, it is obvious that (41) results from (28) by just replacing G 
with | T |. A comparison between the expressions of G and T, that are given by Eqs (29) 
and (42) respectively, shows that their leading terms, for p = 1 and k = 1 respectively, are 
equal. As regards the whole sums we can draw some conclusions in the case of small ip. 
Using the asymptotic formula 13.5.10, U(a,b,ip) = r ^~^ b ) + ^(1 ^ l 1- ^); °f 



Abramowitz 



& Stegun (1970) that holds for small ip and < b < 1 , for a = p — 1 and 6 = 1/2 and 



substituting in (29) we see that G is close to T. 
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3.1. Mass functions and N-body simulations 

On the other hand, large numerical simulations give very important information about 
mass functions. Such simulations, starting from cosmological initial conditions, can find at 
different values of the redshift z, the number density of haloes of given mass M, denoted 
by N(M, z), and the variance of mass at scale M, denoted by S(M). It is well known that 
these quantities are related by the equation: 

N(M, z)dM = -Z&nS, z) ^IdM (43) 

that can be written in the form: 

N(M, z)dM = -^M [2S T(S, z)\ ^^ dM (44) 

We note here that some quantities in the above Eq. can be evaluated from the results of N- 
body simulations. Sheth & Tormen ( 1999 ) showed that the combination M N p^^ d ln [ cr d ^ M z ^ 



has an almost universal behavior , that is independent of redshift and cosmology, a result 



that was confirmed by large numerical simulations as those of Tinker et al. (2008). We 
recall that cr 2 (M, z) is the variance at mass scale M at redshift z. In the linear regime of 
the evolution it obeys the relation a(M,z) = a(M,z = O)D(z) = l.Q86a(M,z = 0)B sc (z). 



Introducing the variable v = B sc (z)y/S, for constant z, we can write 

^ = 2- (45) 
au v 

Assuming a distribution function T v of the variable v and combining, for constant z, the 
fundamental law of probabilities, 

F(S,z)dS = -F„(u)du (46) 

with (45), we get: 

uT v (v) = 2SF{S,z) (47) 



13 




Fig. 1. — Multiplicity functions, vT v {y) satisfying (47). Squares are the measurements from 



N-body simulations of Tinker et al. (2008) and the solid line is the analytical fit to these 
predictions given by the above authors. The line of the analytical fit was derived using (51) 



and A = 178 in Eqs (Bl), (B2), (B3) and (B4) of the appendix of |Tinker et aL| ( |2008[ ), 
for the evaluation of constants A,a,b,c in Eq (51). Dashed-dot-dot line is the prediction 



for the SC model where TiS^z) is given by Eq.37. Dots are the prediction of the Sheth & 



Tormen (2002) approximation for the EC model (Eq. 38). Thin dashes are the prediction 



of the path integral approximation used in this work (Eq.28). In Eqs (28) and (38) we used 
for (5*0,50) = (0,0) . Note that the Eqs (28), (38) and the numerical solution are almost 
indistinguishable from one another.Thick dashes show the numerical solution of (A5), for 



the ellipsoidal barrier, Zhang & Hui (2006). 
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and using (44) we have: 



M 



■N(M, z) 



dM 



(4f 



Pb [z) ■ ■dln(5'- 1 /2 ; 

Taking into account that Bsc(z) and y/S evolve with time in the same way according 
to the linear theory, the quantities v and vT v {y) are time independent. vT v {y) is 
called the "multiplicity function". For SC model, using (37) it is easy to show that 
2SJ z SC {S 1 z) = ^z/exp[-^] while for the EC model of 
write: 



Sheth & Tormen 



(2002) we can 



where 



2SF EC {S,z) 



9{l) 



V2 w 2 [l + g{j)] i 
v-=\wi + =- exp 

'71 V 1 ^ 



1 



w 2 l2 



7 + 



7(7 - 1) 



7(7-1) •••(7-4) 



(49) 



(50) 



2! 5! 
Numerical experiments (see for example Tinker et al. ( 2008[ )) show that vT v {y) is indeed 
nearly a constant. From their N-body simulations the above authors have shown that: 



2SJ-(S, z) n _i, oc iy — A 



S 



+ 1 



e s 



(51) 



where the constants A, a, b and c are given by the Eqs (Bl), (B2), (B3) and (B4) of Tinker 



et al. (2008) 



In what follows, we have assumed a flat model for the Universe with present day 
density parameters f2 m0 = 0.3 and = A/3Hq = 0.7 where A is the cosmological 
constant and H is the present day value of Hubble's constant. We have used the value 
Hq = 100 hKMs _1 Mpc~ 1 and a system of units with m un i t = 1O 12 M /i~ 1 , r un i t = l/i^Mpc 
and a gravitational constant G — 1. At this system of units H /H unit = 1.5276. 
Regarding the power spectrum, we employed the AC DM formula proposed by |Smith et al. 
I (1998). The power spectrum is normalized for <j$ = cr(R = 8/i -1 Mpc) = 0.9. 

In Fig. 1 we have plotted multiplicity functions, vTyiv^. Squares are the predictions from 



N-body simulations of Tinker et al. (2008) and the solid line is the analytical fit to these 
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10 l 



r 



10- 1 H 



CO 



\ 1 0" h 

CO 



10"' - 



V 



Constrained distribution for: 

S =6.8(M =10 12 M so , ar /h) 

z =0.0 

z=0.05 

Dashed line:Path integral approach 
Dotted line: ST02 approximation 
Open squares: Numerical solution 



□ □ 



J L 



10 



15 



20 



Fig. 2. — Constrained distributions, ^(S, z/S , z ) for S = 6.8 that corresponds to M = 
1O 12 /i _1 M , z = 0.05 and z = 0. We used values of S in the range [S , S max ] where 
S m ax = S(10~ 2 M ) ~ 22. All lines are predictions for the ellipsoidal barrier. Dotted line 
corresponds to the ST02 approximation, dashed line to the path approach and squares are 
from the numerical solution of (A5). 
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4.5 




z =0 
z=0.05 



J I I I I L 



25 



Fig. 3.— The relative difference (AT/F) model = {T model - 7 num ) / \Fnum) , where the 
index model stands for the path or the ST02 approach and the index num is for the 
numerical solution, as a function of S. Solid lines stand for the path approach while 
dotted lines for the ST02 approximation. All curves correspond to z — 0.05 and 
z = but are predicted for different ranges of mass. The following values of M, 
lO lo M //i, W n M & /h, 10 12 M G /h, W 13 M G /h, 10 M M G /h correspond to S = 1.18, 3.09, 6.8, 13 
and 22, respectively 
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Fig. 4.— The values of ±J^ c (So, z /S, z) \ Zo as a function of £ (the ratio of the mass of 
the progenitor to the mass of the descendant halo), for S = 6.8. This Fig. is for z = 
and Zq = 3. Dashed lines correspond to the path approach model, dotted lines to the ST02 
approximation while solid lines correspond to the SC model. Upper lines correspond to 
z = and lower lines to z = 3. 
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predictions given by the above authors. Dashed-dot-dot line is the prediction for the SC 



model, where ^(S, z) is given by Eq.37. Dots are the prediction of the Sheth & Tormen 



(2002) approximation for the EC model (Eq. 38). Thin dashes are the prediction of the path 
integral approximation used in this work (Eq.28). Thick dashes show the numerical solution 



of (A5) for the ellipsoidal barrier, Zhang & Hui (2006). We note that the results of Sheth 



& Tormen (2002) approximation and the results of the path integral approach approximate 



the results of N-body simulations better than the results of SC model, especially for values 
of v < 1. This happens not only for mass functions (Sheth & Tormen ( 2002[ )) but for other 



characteristics too as for example the formation times of dark matter haloes (eg. Lin et al. 



( |2003l ), [Hiotelis fc Del Popolol p006| )). 

In Fig. 2 we present constrained distributions for different halo mass and different redshifts. 
These distributions are given by (40), (41) and (28) for the SC, ST02 and the path 
approach models, respectively. Constrained distributions are essential when dealing with 
the formation of dark matter haloes. Important characteristic of dark matter haloes, as 
formation times, merger rates etc., depend on the form of the constrained distributions. 
In this figure -for fixed values of So, z and z ( for z > zo)- the quantity ^(S, z/Sq, Zq) is 
plotted as a function of S (for S > So). This quantity is calculated from both analytical 
approximations and numerical solutions of (A5). It is shown that for values of S close to 
So -that correspond to masses M close to Mo- analytical and numerical methods give the 
same results, but for large values of S - that correspond to values of mass M much smaller 
than mass M - the numerical solutions - assumed to be the exact solutions- give values 
for T that are significantly smaller than those of analytical approximations. In order to 
give a more accurate comparison between the results we calculated the relative difference 



(AJ 7 /J 7 ) mode/ = 



model 



Fnum) I i^num) , where the index model stands for the path or the 



ST02 approach and the index num is for the numerical solution. Both T num and Fmodei 
are constrained distributions (S, z/Sq, z q ). We present results in Fig.3. Solid lines stand for 
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the path approach while dotted lines for the ST02 approximation. All curves correspond to 
z = 0.05 and z = but are predicted for different ranges of mass. The following values 
of M, M = lO lo M //i,M = W u M e /h,M = 1O 12 M //i,M = 10 13 M G /h,M = 10 14 M G /h 
correspond to S — 1.18,3.09,6.8,13 and 22 respectively. The results for other redshifts, 
in the range range [0.3] that we studied, are similar and show that the predictions of the 
path approach are closer to the - exact- curve ( predicted by numerical solutions) than the 
predictions of ST02 approximation. 

We note that the relative difference is an increasing function of S starting from zero at 
S = S . Thus the quantity F x that is defined by the relation F x = J s x J P (S,z/S ,zo)dS, 
gives the fraction of walks that start from the point (So, B(Sq, Zo)) and pass from the point 
(S,B(S,z)) with S in the range [S^S^]. We define as S x the value of S that satisfies 
{AT / JF) mo dei{S x ) = x/100 and thus F x equals to the fraction of walks that agree with 
the exact solution, better than x percent. For all values of x we have S x>pa th > S Xj st and 
F x ^ P ath > F X; st- We have also calculated M x = M(S X ). In the following tables we give some 
characteristic results for x = 10. 

It is clear from the above two tables that path integral approach agrees to the exact results 

Table 1: Results for z = 0, z = 0.05 

Mass S WtS T Sw, P ath Fio,st F 10tPath M WtS T M WtPath 

1O 12 /i _1 M 8.463 9.647 0.9854 0.9894 0.48 x 1O 12 /i _1 M 0.31 x 1O 12 /i _1 M 
W z h~ l M P) 3.970 4.486 0.9623 0.9670 0.51 x l^h~ l M P) 0.35 x lO 13 /*"^ 



U.3IU I.IOU U.»UiO U.3UIU U.dl /\ 1U II IV1 Q V-OO /\ -LU II 1V1 Q 

1.686 1.976 0.9498 0.9570 0.44 x 1O 14 /i" 1 M 0.31 x 1O 14 /i" 1 M 



for a significant larger range of S than ST approximation. Sio, pa th can be by 26 percent 
larger than S 10j st as for example it can be seen in the last row of Table 2. Additionally, 
the results of path approach extended to larger masses. As it can be seen the interval of 
masses [M, Mo] = [M(S), M(Sq)], (note that S > So), with an accuracy better than 10 
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Table 2: Results for zo = 3, z = 3.05 

Mass S Wt sr Sio )Pa th F Wi st FiQ^p a th M Wt sT M WtPath 

W^h^MQ 9.025 10.536 0.9691 0.9754 0.38 x 10 12 h^M & 0.22 x W 2 h~ l M Q 

IO^/i-^q 4.458 5.383 0.9532 0.9613 0.36 x IO^/^Mq 0.21 x IO^H^Mq 

10 14 h- 1 M Q 1.954 2.475 0.9392 0.9504 0.31 x 1O 14 /i~ 1 M 0.18 x 1O 14 /i~ 1 M 

percent becomes larger for the path approach since it extends to smaller values of mass. On 
the other hand, the same Tables show that the improvement in F 10 is not significant. This 
improvement seems to be an increasing function of M and is less than 1.2 percent for the 
range of parameters we studied. Similar results were found for other values of x too. 



4. Merger rates 



Using Bayes rule we write: 

F c (S ,z /S,z) 
Using (29), (41) or (42) we can write: 

T c (S ,z /S,z) = 



7XS,z/S ,zo).F(,So,zo) 
T(S,z) 



3 / 2 jS B(S,z)-S5 ] 2 



where 



$ m (S, So, z, 5 ) = 



So(S-S ) 

G m (S, Sq, z, <5o)G m (So, 0, zq, 0) 



(52) 



e 2ss 0( ss 0) $ m (5,So,z,<5 ) (53) 



(54) 



G m (S, 0,^,0) 

The index m states for the following models: SC model, ST-EC model and the path model, 
P-EC. Thus: 

G SC (S, S , z, 5 ) = B sc (z) -5 = B sc (z) - B sc (z ) (55) 



and 



Gst-ec{S, S , z, S Q ) =| T(S, zj S , z Q ) 



(56) 
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Fig. 5. — Mean merger rates for M = 10 11 M Q /h and for three different redshifts z = 0, z = 
1 and z = 3. Dashed lines (from top to bottom) correspond to the path approach for the 
above redshifts, respectively. Solid lines are the predictions of ST02 model and dotted lines 
are the predictions of SC model. 
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while G P _ EC (S, S , z, S ) is given by (29). Differentiating (53) with respect to z we have for 
the various models the following relations: 
For the SC model: 



dz 



FscWoi z o/ S, z) 



>0 



2tt 
So 



1 3/2 



B 2 sc {z) 



Sq(S — So 
1 - 



cxp 



[S B sc (z) - S5 ] 
2S Sq(S — iSq) 



x 



S \ S B S c(z) - SS 



Bsc(z) 



S(S - So) 



dB S c(z) 
dz 



(57) 



For the ST - EC model, 



&FsT-Ec{ S 0> z o/ S > z ) _ 7 /q 7 \v 
-j- — J- ST-ECWO, z ol O, Z) X 

SB(S , z ) - S B(S, z) dB(S, z) T(S, z) d (T(S, z/S , z ) 



S(S - S ) dz T(S, z/S , zo) dz \ T(S, z 

and for the P — EC model: 

dJ r p_ EC (So, zoj S, z) 



(58) 



dz 



Fp-EciSo, Zo/S, z) x 



SB(So,z ) - S B(S,z) dB(S,z) G P . EC (S, 0, z, 0) d fG P _ EC (S,S ,z,S ) 



S(S - S ) dz Gp- EC (S, S 0: z, S ) dz \ G P - EC (S, 0, z, 0) 



(59) 



In Fig 4 we present -^J^sci^O: z o/S,z) \ zo as a function of ^ for S = 6.8 ( that corresponds 



dz' 

U2 i 



to Mo = 10 Mq/K) and zo = 0, zo = 3. Details about the particular values of the 
parameters are written in the figure caption. The results show some interesting features: 

a) SC model gives completely different results. Only for large values of £, that means 
for large progenitors, the results of SC model are somehow close to those predicted by 
the models that use the ellipsoidal barrier. For small values of £, that means for small 
progenitors SC and EC models lead to completely different results. 

b) From the study of the same quantity for different values of the parameters mass M and 
redshift zo we conclude that the predictions of path approach and the ST02 approximation 
are close for any redshift and mass in the range [0,3.] and [10 12 M o //i, 1O 14 M //i] 
respectively. 
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We define as merger rate at redshift z the quantity R, 



df dJ 70 
R(M M /z )dM = J-(M , zq/M, z) \ z=zo dM = — (S , z /S, z) \ z=zo dS (60) 

R(M — > M /z ) equals to a fraction of the (total) mass that belongs to haloes of masses M. 
It is this fraction of mass that merges instantaneously to form at z haloes of mass M . The 
quantity R(M — > M /z )f(M,z )dM expresses the mass that belongs to haloes of mass 
M, M + dM, and merges instantaneously to form at z haloes of mass M , as a fraction of 
the total mass of the Universe. Multiplying R(M M /z )f(M, z )dM by V un p b /M and 
dividing by V un f(M , zo)pb/M where V un is the volume of the universe, we find: 

N m (M , z )dz M /(M , z ) 
In the above relation N m (M, zq) is the number of haloes of mass in the range M, M + dM 
that merge at z and form haloes of mass M while N m (M , z ) is the number of haloes of 
mass M present at z . Using (33), we write: 

N m (M,Zo) M,7MdSd 

iv m (A/ 0>2o) d, = m ns^y dis J|4 2 ° /s ' 2)dM (62) 

In a merger procedure as above, the halo resulting from mergers between smaller ones, a 
halo with mass M in our notation, is called a descendant halo. Haloes with smaller masses 
that merge to form M are called progenitors. It is convenient for comparing merger rates 
for different descendant haloes to use, instead of the masses M and M , the ratio £ = M/M 
and thus (61) is written: 



B (rn) = 



iV m (M,z ) M 2 ^(5,^o) 



dJ- c (5o,^/g,^) 
dz 



,lS (63) 



2 dM 



M °'*° vw iV m (M , ^ )d^d£ M F(S ,z ) 
and it expresses the mean merger rate. It measures the mean number of mergers per halo, 
per unit redshift, for descendant haloes of mass Mo with mass ratio £ as defined above. 
The mean merger rates for the three models are: 

'dBsc(z) 



dM,SC/M _ 1 M 0(C q \-3/2i dS I 

B *»*> {0 ~ 7^Tf {S ~ So) 'dM 1 



dz 



(64) 

20 
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Fig. 6.— As in Fig.5 but for M = 1O 14 M //i. 
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B 



SD M< i ,z Q \S>) 

(m),P-EC , (.s 
M ,z 



M fQ /C^3/2| T(S,Z ) 



T'('S'o) zo) 



M 2 



SV>\ 3 ^ 2 Gp-ec(S, 0, 0) AB 



ST-EC 

dz 



d T c 

U>A P-EC 



dS 
dM 



M \ S J Gp- EC (S o ,0,z o ,0) 



dz 



2() 



dS 
dM 



(65) 
(66) 



where AB = — ^ 



B| c (So,*o) 



fl| (S,*D) 
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In Figs 5 and 6 we present mean merger rates derived by Eqs 64, 65 and 66 for different 
redshifts and masses of descendant haloes. Details are written in figures captions. It is 
clear that for small haloes the predicted merger rates of SC model are far from those of 
the EC model that predicted by path approach (see Fig. 5). For increasing halo masses 
the agreement becomes better. We note that analyzing the results of N-body simulations 



Fakhouri & Ma (2008) and Stewart et al. (2009) found analytical approximations for 



merger rates of dark matter haloes but their definition of £ differs to the one used above. In 
their papers, £ is not defined as £ = M/M but as £ = M/Mmmp, where Mmmp is the most 
massive progenitor during the merge process. With this definition of £ merger rates cannot 
be found analytically. Their derivation requires the construction of merger-trees, that is 



a complex process with its own difficulties, (e.g. Neinstein & Dekel (2008)). In Hiotelis 



(2011), using merger-tress, it was shown that SC model approximates better merger rates 
for small haloes while merger rates of larger haloes are approximated better from models 
that use the ellipsoidal barrier (see Fig. 9 therein). 

A direct detailed comparison of mean merger rates derived by Eqs 64, 65 and 66 -where 
£ = M/Mq- with those predicted from N-body simulations is the next step in our research. 



5. Conclusions and Discussion 

The power tool of path integrals is used in order to approximate the problem of the 
formation of dark matter haloes. The initial density field is assumed to be Gaussian and 
it is smoothed by a filter that is sharp in k space. Using these assumptions, the results 
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show that decreasing the radius of the filter, the smoothed density evolves according to 
a Langevin equation of a form that describes a Markovian stochastic process. Thus, it 
was able to calculate first crossing distributions for two -well known in the literature- 
namely a flat barrier corresponding to the Spherical Collapse model and a moving barrier 
corresponding to the Ellipsoidal Collapse Model. Mass functions of dark matter haloes, 
predicted by the above first crossing distributions, are compared with numerical solutions 
-that are considered as exact- and with the predictions of N-body simulations. It is shown 
that for the EC model -that is more promising than the SC model- the prediction of path 



approach is close to the results of Sheth & Tormen (2002). 

Although in the calculation of T we used only two terms from an expansion (see Eqs. 28 
and 29), the resulting constrained first crossing distributions from the path approach are 
closer to the ones predicted by the numerical solutions than the predictions of |Sheth fc 



Tormen (2002) are. We note that the infinite sum in (29) is accurately approximated using 



up to 12 terms. This shows, that the analytical formula given in Eq. 28, that resulted from 
the analytical procedure described in the text, improves the empirical formula (see Eqs. 41 



and 42) of Sheth & Tormen (2002). We note that the predictions of the analytical formula 
of ST02 are predicted using about 6-8 terms of the sum. 

Merger rates have also been calculated. Merger rates resulting from the path integral 



approach are close to those predicted by Sheth & Tormen (2002) approximation and no 
significant improvement appears. The use of additional terms in Eqs. 28 and 29 -that 
could bring the results of path approach closer to that of the numerical solutions of the 
integral equation (A5)- as well as a detailed comparison with the merger rates resulting 
from N-body simulations are clearly required. Both actions are in progress. 
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A. Appendix A 



We denote by P(S,5/So,So)dS the probability a trajectory that starts from the point 
(So, do) passes at "time" S between 5, 5 + d<J without crossing the barrier B(S, z) between 
So and S. On the other hand J-(S, z/Sq, 5o)dS = J-"(S, B(S, z)/Sq, So)dS equals to the 
probability a trajectory that passes from the point (So, So) crosses the barrier of height 
B(S, z) > 8q, for the first time, between S, S + dS. Consequently, Pi given by 



Pi 



F(S',B(S',z)/S ,5o)dS' 



(Al) 



So 



is the probability the trajectory has crossed the barrier before S while P2, 

"B(S,z) 



Po 



P(S,6/S ,5 )d5 



(A2) 



is the probability the trajectory was always at values smaller than B(S,z) and thus it has 
not crossed the barrier for "time" < S. Obviously Pi + P2 = 1- We also assume that the 
transition probability Pq, in the absence of any barrier, is a normal Gaussian given by: 



Po(S,S/S Q ,S 



1 



cxp 



VMS - s 

The presence of the barrier amplifies P in the following way: 

-s 



(S - ^o) 2 
2(S - S ) 



(A3) 



P(S, 6 /So, 5 ) = P (S, 5/So, (Jo) - / F[S', B(S', z)/S Q , 5 ]Po[S, 8/S', B(S', z)]dS' (A4) 

J So 



Using (A3) and (A4) it can be proved (Zhang & Hui 2006), that for an arbitrary barrier 



T satisfies the following integral equation: 

T[S, B(S, z)/S , So) = gi (S, 60, So) + f g 2 (S, S')T[S', B(S', z)/S , S ]dS' 

J So 



where: 



gi(S,S ,S ) 



B(S,z)-8o _ 2 dB(S,z) 



92(S, S') 



S-S dS 
dB(S,z) B(S,z) - B(S',z 



dS 



S-S' 



Po[S,B(S,z)/So,So] 
Po[S,B(S,z)/S'B(S',z)} 



(A5) 

(A6) 
(A7) 
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In the case of a linear barrier Eq.(A5) admits an analytic solution. If B(S, z) = p c (z)+q c (z)S, 
where p c and q c are functions of the redshift z, the solution is written: 



JF[S, z/So, d \ = j=. = == exp 



[B(S,z)-5 ] 2 
2(S - S ) 



(A8) 



y/2n(S-S ) 3 

Thus, the spherical model which is of the form B(S,z) = B(z) = 1.686/ D(z) leads to the 
solution: 



Fsc{S,z/So,zq) 



B{z) - 5 



exp 



[B(z) - 5 
2(S - So) 



(A9) 



^2n(S-S ) 3 

Eq. (A5) admits a simple numerical solution. First, we use a grid of points Sk = S + kAS 
for k = 0, 1...N and AS = (S - S )/N (thus S N = S), for the interval [S , S}. Then, using 
the trapezoidal rule, we write (A5) as: 

j=i A<7 

(AlO) 



AS i \ AS 

?i = 9i(Si,S ,S ) + — ^2g 2 (S i ,S j 7r)\Fj + ^-1} 



2 2 

j'=i 

where we set Tk = TiS^, z/S , z ). Solving for J 7 ; we have the solution of the integral 
equation: 



Fi = 9i(S» Sp, Sp) + f Eg!' ftCfl, g; - + Jj-i] 

1 — 92(81, Si — 4r) 
The above Eq. holds for i > 2 while for i = or % = 1 we have J-q — and 

T\ = g 1 (S 1 ,5 ,S )/[l - ^fg 2 (S 1 ,S 1 - ^f)], respectively. 

We also used an iterative method for the solution of (A5). We denote by the I th 
approximation for the value J 7 ,. We use the initial conditions F- ^ = gi(S i7 5 , S ) and the 
iterative method: 



= 9^,60, S ) + ^5>( 5i ,S, - ^)[jf + r/i, 



(I) 



2 2 



(A12) 



The iterations stop when for all values of % — 1, 2..N, the condition: 



(A13) 



We used e — 10 5 . We note that the numerical solution of (A5) resulting from any of 
the above two methods is very sensitive to the number of grid points N. We used a large 
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number of grid points N, N = 10 5 , to divide the interval [S , S max ] where S = S(M ) and 
Smax = S(10~ 2 M ) so as the two numerical methods give identical results. The accuracy of 
these numerical results was checked by comparing them with the exact solutions that are 
available from the SC model. The agreement was found almost exact. 



